/* Copyright (c) 2021, National University of Defense Technology. All rights reserved.
 *
 * Licensed under the Apache License, Version 2.0 (the "License");
 * you may not use this file except in compliance with the License.
 * You may obtain a copy of the License at
 * 
 *     http://www.apache.org/licenses/LICENSE-2.0
 * 
 * Unless required by applicable law or agreed to in writing, software
 * distributed under the License is distributed on an "AS IS" BASIS,
 * WITHOUT WARRANTIES OR CONDITIONS OF ANY KIND, either express or implied.
 * See the License for the specific language governing permissions and
 * limitations under the License.
 */

#include <set>
#include <map>
#include <stdlib.h>
#include "parcsr_coarsen.h"

namespace YHAMG
{

#define COARSE 1
#define FINE -1
#define UNASSIGNED 0

#define MPI_TAG 400

static inline void Swap(int& a, int& b)
{
	int temp = a;
	a = b;
	b = temp;
}

#define RANDOM_MAX 32767

static inline unsigned int Random(unsigned int& x)
{
	x ^= x << 13;
	x ^= x >> 17;
	x ^= x << 5;
	return x & 32767;
}

int ParCSR_HMIS_Coarsening(const ParCSR_Matrix& A, const bool* LocLoc_s, const bool* LocExt_s, int* CF_map)
{
	MPI_Comm comm = A.comm;

	int rank;
	MPI_Comm_rank(comm, &rank);

	int n = A.LocLoc.size[0];
	int Ext_size = A.LocExt.size[1];

	int* LocLoc_rowptr = A.LocLoc.rowptr;
	int* LocLoc_colind = A.LocLoc.colind;
	int* LocExt_rowptr = A.LocExt.rowptr;
	int* LocExt_colind = A.LocExt.colind;

	int num_interface_rows = 0;
	int* interface_rows = new int[n];
	for (int i = 0; i < n; ++i)
		if (LocExt_rowptr[i + 1] > LocExt_rowptr[i])
			interface_rows[num_interface_rows++] = i;

	int num_neighbors = A.num_neighbors;
	int* neighbor_ranks = A.neighbor_ranks;
	int* recv_ptr = A.Ext_ptr;
	int* send_ptr = A.send_ptr;
	int* send_list = A.send_list;
	int* send_buffer = new int[send_ptr[num_neighbors]];

	int* weights = new int[n];
	int* Ext_weights = new int[Ext_size];

	int* Stp = new int[n + 1];
	for (int i = 0; i <= n; ++i)
		Stp[i] = 0;

	for (int j = LocLoc_rowptr[0]; j < LocLoc_rowptr[n]; ++j)
		if (LocLoc_s[j]) ++Stp[LocLoc_colind[j]];

	for (int i = 0; i < n; ++i)
		Stp[i + 1] += Stp[i];

	int* Sti = new int[Stp[n]];

	for (int i = n - 1; i >= 0; --i)
		for (int j = LocLoc_rowptr[i + 1] - 1; j >= LocLoc_rowptr[i]; --j)
			if (LocLoc_s[j]) Sti[--Stp[LocLoc_colind[j]]] = i;

	for (int i = 0; i < n; ++i)
		weights[i] = Stp[i + 1] - Stp[i];

	int num_coarse = 0;
	int num_fine = 0;

	for (int i = 0; i < n; ++i)
	{
		if (weights[i] == 0)
		{
			CF_map[i] = FINE;
			++num_fine;
		}
		else
			CF_map[i] = UNASSIGNED;
	}

	while (num_coarse + num_fine < n)
	{
		for (int i = 0; i < n; ++i)
		{
			int r = i;

			if (CF_map[r] == UNASSIGNED)
			{
			L0:
				for (int j = LocLoc_rowptr[r]; j < LocLoc_rowptr[r + 1]; ++j)
				{
					int jcol = LocLoc_colind[j];

					if (CF_map[jcol] == UNASSIGNED && weights[r] < weights[jcol])
					{
						r = jcol;
						goto L0;
					}
				}

				for (int j = Stp[r]; j < Stp[r + 1]; ++j)
				{
					int jcol = Sti[j];

					if (CF_map[jcol] == UNASSIGNED && weights[r] < weights[jcol])
					{
						r = jcol;
						goto L0;
					}
				}

				CF_map[r] = COARSE;
				++num_coarse;

				for (int j = Stp[r]; j < Stp[r + 1]; ++j)
				{
					int jcol = Sti[j];
					if (CF_map[jcol] == UNASSIGNED)
					{
						CF_map[jcol] = FINE;
						++num_fine;

						for (int k = LocLoc_rowptr[jcol]; k < LocLoc_rowptr[jcol + 1]; ++k)
						{
							int kcol = LocLoc_colind[k];

							if (LocLoc_s[k] && CF_map[kcol] == UNASSIGNED)
								++weights[kcol];
						}
					}
				}

				for (int j = LocLoc_rowptr[r]; j < LocLoc_rowptr[r + 1]; ++j)
				{
					int jcol = LocLoc_colind[j];

					if (LocLoc_s[j] && CF_map[jcol] == UNASSIGNED)
					{
						if (--weights[jcol] == 0)
						{
							CF_map[jcol] = FINE;
							++num_fine;

							for (int k = LocLoc_rowptr[jcol]; k < LocLoc_rowptr[jcol + 1]; ++k)
							{
								int kcol = LocLoc_colind[k];

								if (LocLoc_s[k] && CF_map[kcol] == UNASSIGNED)
									++weights[kcol];
							}
						}
					}
				}
			}
		}
	}

	int lambda = RANDOM_MAX + 1;
	static unsigned int x = rank + 2671;

	for (int i = 0; i < num_interface_rows; ++i)
	{
		int irow = interface_rows[i];
		weights[irow] = (Stp[irow + 1] - Stp[irow]) * lambda + Random(x);
	}

	for (int i = 0; i < Ext_size; ++i)
		Ext_weights[i] = 0;

	for (int j = LocExt_rowptr[0]; j < LocExt_rowptr[n]; ++j)
		if (LocExt_s[j]) Ext_weights[LocExt_colind[j]] += lambda;

	MPI_Request* recv_requests = new MPI_Request[num_neighbors];
	MPI_Request* send_requests = new MPI_Request[num_neighbors];
	MPI_Status status;

	for (int r = 0; r < num_neighbors; ++r)
		if (send_ptr[r + 1] > send_ptr[r])
			MPI_Irecv(send_buffer + send_ptr[r], send_ptr[r + 1] - send_ptr[r], MPI_INT, neighbor_ranks[r], MPI_TAG, comm, recv_requests + r);

	for (int r = 0; r < num_neighbors; ++r)
		if (recv_ptr[r + 1] > recv_ptr[r])
			MPI_Isend(Ext_weights + recv_ptr[r], recv_ptr[r + 1] - recv_ptr[r], MPI_INT, neighbor_ranks[r], MPI_TAG, comm, send_requests + r);

	for (int r = 0; r < num_neighbors; ++r)
	{
		if (send_ptr[r + 1] > send_ptr[r])
		{
			MPI_Wait(recv_requests + r, &status);
			for (int i = send_ptr[r]; i < send_ptr[r + 1]; ++i)
				weights[send_list[i]] += send_buffer[i];
		}
	}

	for (int r = 0; r < num_neighbors; ++r)
		if (recv_ptr[r + 1] > recv_ptr[r])
			MPI_Wait(send_requests + r, &status);

	for (int r = 0; r < num_neighbors; ++r)
		if (recv_ptr[r + 1] > recv_ptr[r])
			MPI_Irecv(Ext_weights + recv_ptr[r], recv_ptr[r + 1] - recv_ptr[r], MPI_INT, neighbor_ranks[r], MPI_TAG, comm, recv_requests + r);

	for (int r = 0; r < num_neighbors; ++r)
	{
		if (send_ptr[r + 1] > send_ptr[r])
		{
			for (int i = send_ptr[r]; i < send_ptr[r + 1]; ++i)
				send_buffer[i] = weights[send_list[i]];
			MPI_Isend(send_buffer + send_ptr[r], send_ptr[r + 1] - send_ptr[r], MPI_INT, neighbor_ranks[r], MPI_TAG, comm, send_requests + r);
		}
	}

	for (int r = 0; r < num_neighbors; ++r)
		if (recv_ptr[r + 1] > recv_ptr[r])
			MPI_Wait(recv_requests + r, &status);

	for (int r = 0; r < num_neighbors; ++r)
		if (send_ptr[r + 1] > send_ptr[r])
			MPI_Wait(send_requests + r, &status);

	delete[] Stp;
	delete[] Sti;

	for (int i = 0; i < num_interface_rows; ++i)
	{
		int irow = interface_rows[i];
		--(CF_map[irow] == COARSE ? num_coarse : num_fine);
		CF_map[irow] = UNASSIGNED;
	}

	for (int i = 0; i < num_interface_rows; ++i)
	{
		int irow = interface_rows[i];

		if (weights[irow] < lambda)
		{
			CF_map[irow] = FINE;
			++num_fine;
			continue;
		}

		for (int j = LocLoc_rowptr[irow]; j < LocLoc_rowptr[irow + 1]; ++j)
		{
			int jcol = LocLoc_colind[j];

			if (LocLoc_s[j] && CF_map[jcol] == COARSE)
			{
				CF_map[irow] = FINE;
				++num_fine;
				break;
			}
		}

	}

	int* Ext_CF_map = new int[Ext_size];

	for (int r = 0; r < num_neighbors; ++r)
		if (recv_ptr[r + 1] > recv_ptr[r])
			MPI_Irecv(Ext_CF_map + recv_ptr[r], recv_ptr[r + 1] - recv_ptr[r], MPI_INT, neighbor_ranks[r], MPI_TAG, comm, recv_requests + r);

	for (int r = 0; r < num_neighbors; ++r)
	{
		if (send_ptr[r + 1] > send_ptr[r])
		{
			for (int i = send_ptr[r]; i < send_ptr[r + 1]; ++i)
				send_buffer[i] = CF_map[send_list[i]];
			MPI_Isend(send_buffer + send_ptr[r], send_ptr[r + 1] - send_ptr[r], MPI_INT, neighbor_ranks[r], MPI_TAG, comm, send_requests + r);
		}
	}

	for (int r = 0; r < num_neighbors; ++r)
		if (recv_ptr[r + 1] > recv_ptr[r])
			MPI_Wait(recv_requests + r, &status);

	for (int r = 0; r < num_neighbors; ++r)
		if (send_ptr[r + 1] > send_ptr[r])
			MPI_Wait(send_requests + r, &status);

	int* neighbor_unassigned = new int[num_neighbors];
	for (int r = 0; r < num_neighbors; ++r)
		neighbor_unassigned[r] = 1;

	int LESS = 0;
	int GREATER = 2;

	while (1)
	{
		int num_unassigned = n - num_coarse - num_fine;

		for (int r = 0; r < num_neighbors; ++r)
			if (neighbor_unassigned[r] != -1)
				MPI_Irecv(neighbor_unassigned + r, 1, MPI_INT, neighbor_ranks[r], MPI_TAG, comm, recv_requests + r);

		for (int r = 0; r < num_neighbors; ++r)
			if (neighbor_unassigned[r] != -1)
				MPI_Send(&num_unassigned, 1, MPI_INT, neighbor_ranks[r], MPI_TAG, comm);

		for (int r = 0; r < num_neighbors; ++r)
		{
			if (neighbor_unassigned[r] != -1)
			{
				MPI_Wait(recv_requests + r, &status);
				if (neighbor_unassigned[r] == 0)
					neighbor_unassigned[r] = -1;
			}
		}

		if (num_unassigned == 0) break;

		Swap(LESS, GREATER);

		for (int i = 0; i < num_interface_rows; ++i)
		{
			int irow = interface_rows[i];

			if (CF_map[irow] == GREATER)
			{
				for (int j = LocLoc_rowptr[irow]; j < LocLoc_rowptr[irow + 1]; ++j)
				{
					int jcol = LocLoc_colind[j];

					if (LocLoc_s[j] && CF_map[jcol] == GREATER)
					{
						if (weights[jcol] > weights[irow])
						{
							CF_map[irow] = LESS;
							break;
						}
						else
							CF_map[jcol] = LESS;
					}
				}
			}
		}

		for (int i = 0; i < num_interface_rows; ++i)
		{
			int irow = interface_rows[i];

			if (CF_map[irow] == GREATER)
			{
				for (int j = LocExt_rowptr[irow]; j < LocExt_rowptr[irow + 1]; ++j)
				{
					int jcol = LocExt_colind[j];

					if (LocExt_s[j] && Ext_CF_map[jcol] == GREATER)
					{
						if (Ext_weights[jcol] > weights[irow])
						{
							CF_map[irow] = LESS;
							break;
						}
					}

				}
			}
		}

		for (int r = 0; r < num_neighbors; ++r)
			if (recv_ptr[r + 1] > recv_ptr[r] && neighbor_unassigned[r] != -1)
				MPI_Irecv(Ext_CF_map + recv_ptr[r], recv_ptr[r + 1] - recv_ptr[r], MPI_INT, neighbor_ranks[r], MPI_TAG, comm, recv_requests + r);

		for (int r = 0; r < num_neighbors; ++r)
		{
			if (send_ptr[r + 1] > send_ptr[r] && neighbor_unassigned[r] != -1)
			{
				for (int i = send_ptr[r]; i < send_ptr[r + 1]; ++i)
					send_buffer[i] = CF_map[send_list[i]];
				MPI_Isend(send_buffer + send_ptr[r], send_ptr[r + 1] - send_ptr[r], MPI_INT, neighbor_ranks[r], MPI_TAG, comm, send_requests + r);
			}
		}

		for (int r = 0; r < num_neighbors; ++r)
			if (recv_ptr[r + 1] > recv_ptr[r] && neighbor_unassigned[r] != -1)
				MPI_Wait(recv_requests + r, &status);

		for (int r = 0; r < num_neighbors; ++r)
			if (send_ptr[r + 1] > send_ptr[r] && neighbor_unassigned[r] != -1)
				MPI_Wait(send_requests + r, &status);

		for (int i = 0; i < num_interface_rows; ++i)
		{
			int irow = interface_rows[i];

			if (CF_map[irow] == GREATER)
			{
				CF_map[irow] = COARSE;
				++num_coarse;
			}
		}

		for (int i = 0; i < num_interface_rows; ++i)
		{
			int irow = interface_rows[i];

			if (CF_map[irow] == LESS)
			{
				for (int j = LocLoc_rowptr[irow]; j < LocLoc_rowptr[irow + 1]; ++j)
				{
					if (LocLoc_s[j] && CF_map[LocLoc_colind[j]] == COARSE)
					{
						CF_map[irow] = FINE;
						++num_fine;
						break;
					}
				}
			}
		}

		for (int i = 0; i < num_interface_rows; ++i)
		{
			int irow = interface_rows[i];

			if (CF_map[irow] == LESS)
			{
				for (int j = LocExt_rowptr[irow]; j < LocExt_rowptr[irow + 1]; ++j)
				{
					if (LocExt_s[j] && Ext_CF_map[LocExt_colind[j]] == GREATER)
					{
						CF_map[irow] = FINE;
						++num_fine;
						break;
					}
				}
			}
		}

		for (int r = 0; r < num_neighbors; ++r)
			if (recv_ptr[r + 1] > recv_ptr[r] && neighbor_unassigned[r] != -1)
				MPI_Irecv(Ext_CF_map + recv_ptr[r], recv_ptr[r + 1] - recv_ptr[r], MPI_INT, neighbor_ranks[r], MPI_TAG, comm, recv_requests + r);

		for (int r = 0; r < num_neighbors; ++r)
		{
			if (send_ptr[r + 1] > send_ptr[r] && neighbor_unassigned[r] != -1)
			{
				for (int i = send_ptr[r]; i < send_ptr[r + 1]; ++i)
					send_buffer[i] = CF_map[send_list[i]];
				MPI_Isend(send_buffer + send_ptr[r], send_ptr[r + 1] - send_ptr[r], MPI_INT, neighbor_ranks[r], MPI_TAG, comm, send_requests + r);
			}
		}

		for (int r = 0; r < num_neighbors; ++r)
			if (recv_ptr[r + 1] > recv_ptr[r] && neighbor_unassigned[r] != -1)
				MPI_Wait(recv_requests + r, &status);

		for (int r = 0; r < num_neighbors; ++r)
			if (send_ptr[r + 1] > send_ptr[r] && neighbor_unassigned[r] != -1)
				MPI_Wait(send_requests + r, &status);
	}

	delete[] recv_requests;
	delete[] send_requests;
	delete[] send_buffer;
	delete[] weights;
	delete[] interface_rows;
	delete[] Ext_weights;
	delete[] Ext_CF_map;
	delete[] neighbor_unassigned;

	return num_coarse;
}

int ParCSR_PMIS_Coarsening(const ParCSR_Matrix& A, const bool* LocLoc_s, const bool* LocExt_s, int* CF_map)
{
	MPI_Comm comm = A.comm;

	int rank;
	MPI_Comm_rank(comm, &rank);

	int n = A.LocLoc.size[0];
	int Ext_size = A.LocExt.size[1];

	int* LocLoc_rowptr = A.LocLoc.rowptr;
	int* LocLoc_colind = A.LocLoc.colind;
	int* LocExt_rowptr = A.LocExt.rowptr;
	int* LocExt_colind = A.LocExt.colind;

	int num_interface_rows = 0;
	int* interface_rows = new int[n];
	for (int i = 0; i < n; ++i)
		if (LocExt_rowptr[i + 1] > LocExt_rowptr[i])
			interface_rows[num_interface_rows++] = i;

	int num_neighbors = A.num_neighbors;
	int* neighbor_ranks = A.neighbor_ranks;
	int* recv_ptr = A.Ext_ptr;
	int* send_ptr = A.send_ptr;
	int* send_list = A.send_list;
	int* send_buffer = new int[send_ptr[num_neighbors]];

	int* weights = new int[n];
	int* Ext_weights = new int[Ext_size];

	MPI_Request* recv_requests = new MPI_Request[num_neighbors];
	MPI_Request* send_requests = new MPI_Request[num_neighbors];
	MPI_Status status;

	int lambda = RANDOM_MAX + 1;
	static unsigned int x = rank + 2671;

	for (int i = 0; i < n; ++i)
		weights[i] = Random(x);

	for (int i = 0; i < Ext_size; ++i)
		Ext_weights[i] = 0;

	for (int j = LocLoc_rowptr[0]; j < LocLoc_rowptr[n]; ++j)
		if (LocLoc_s[j]) weights[LocLoc_colind[j]] += lambda;

	for (int j = LocExt_rowptr[0]; j < LocExt_rowptr[n]; ++j)
		if (LocExt_s[j]) Ext_weights[LocExt_colind[j]] += lambda;

	for (int r = 0; r < num_neighbors; ++r)
		if (send_ptr[r + 1] > send_ptr[r])
			MPI_Irecv(send_buffer + send_ptr[r], send_ptr[r + 1] - send_ptr[r], MPI_INT, neighbor_ranks[r], MPI_TAG, comm, recv_requests + r);

	for (int r = 0; r < num_neighbors; ++r)
		if (recv_ptr[r + 1] > recv_ptr[r])
			MPI_Isend(Ext_weights + recv_ptr[r], recv_ptr[r + 1] - recv_ptr[r], MPI_INT, neighbor_ranks[r], MPI_TAG, comm, send_requests + r);

	for (int r = 0; r < num_neighbors; ++r)
	{
		if (send_ptr[r + 1] > send_ptr[r])
		{
			MPI_Wait(recv_requests + r, &status);
			for (int i = send_ptr[r]; i < send_ptr[r + 1]; ++i)
				weights[send_list[i]] += send_buffer[i];
		}
	}

	for (int r = 0; r < num_neighbors; ++r)
		if (recv_ptr[r + 1] > recv_ptr[r])
			MPI_Wait(send_requests + r, &status);

	for (int r = 0; r < num_neighbors; ++r)
		if (recv_ptr[r + 1] > recv_ptr[r])
			MPI_Irecv(Ext_weights + recv_ptr[r], recv_ptr[r + 1] - recv_ptr[r], MPI_INT, neighbor_ranks[r], MPI_TAG, comm, recv_requests + r);

	for (int r = 0; r < num_neighbors; ++r)
	{
		if (send_ptr[r + 1] > send_ptr[r])
		{
			for (int i = send_ptr[r]; i < send_ptr[r + 1]; ++i)
				send_buffer[i] = weights[send_list[i]];
			MPI_Isend(send_buffer + send_ptr[r], send_ptr[r + 1] - send_ptr[r], MPI_INT, neighbor_ranks[r], MPI_TAG, comm, send_requests + r);
		}
	}

	for (int r = 0; r < num_neighbors; ++r)
		if (recv_ptr[r + 1] > recv_ptr[r])
			MPI_Wait(recv_requests + r, &status);

	for (int r = 0; r < num_neighbors; ++r)
		if (send_ptr[r + 1] > send_ptr[r])
			MPI_Wait(send_requests + r, &status);

	int num_coarse = 0;
	int num_fine = 0;

#ifdef USE_OPENMP
#pragma omp parallel for reduction(+:num_fine) schedule(guided)
#endif
	for (int i = 0; i < n; ++i)
	{
		if (weights[i] < lambda)
		{
			CF_map[i] = FINE;
			++num_fine;
		}
		else
			CF_map[i] = UNASSIGNED;
	}

	int* Ext_CF_map = new int[Ext_size];
	for (int i = 0; i < Ext_size; ++i)
		Ext_CF_map[i] = Ext_weights[i] < lambda ? FINE : UNASSIGNED;

	int* neighbor_unassigned = new int[num_neighbors];
	for (int r = 0; r < num_neighbors; ++r)
		neighbor_unassigned[r] = 1;

	int LESS = 0;
	int GREATER = 2;

	while (1)
	{
		int num_unassigned = n - num_coarse - num_fine;

		for (int r = 0; r < num_neighbors; ++r)
			if (neighbor_unassigned[r] != -1)
				MPI_Irecv(neighbor_unassigned + r, 1, MPI_INT, neighbor_ranks[r], MPI_TAG, comm, recv_requests + r);

		for (int r = 0; r < num_neighbors; ++r)
			if (neighbor_unassigned[r] != -1)
				MPI_Send(&num_unassigned, 1, MPI_INT, neighbor_ranks[r], MPI_TAG, comm);

		for (int r = 0; r < num_neighbors; ++r)
		{
			if (neighbor_unassigned[r] != -1)
			{
				MPI_Wait(recv_requests + r, &status);
				if (neighbor_unassigned[r] == 0)
					neighbor_unassigned[r] = -1;
			}
		}

		if (num_unassigned == 0) break;

		Swap(LESS, GREATER);

#ifdef USE_OPENMP
#pragma omp parallel for schedule(guided)
#endif
		for (int i = 0; i < n; ++i)
		{
			if (CF_map[i] == GREATER)
			{
				for (int j = LocLoc_rowptr[i]; j < LocLoc_rowptr[i + 1]; ++j)
				{
					int jcol = LocLoc_colind[j];

					if (LocLoc_s[j] && CF_map[jcol] == GREATER)
					{
						if (weights[jcol] > weights[i])
						{
							CF_map[i] = LESS;
							break;
						}
						else
							CF_map[jcol] = LESS;
					}
				}
			}
		}

#ifdef USE_OPENMP
#pragma omp parallel for schedule(guided)
#endif
		for (int i = 0; i < num_interface_rows; ++i)
		{
			int irow = interface_rows[i];

			if (CF_map[irow] == GREATER)
			{
				for (int j = LocExt_rowptr[irow]; j < LocExt_rowptr[irow + 1]; ++j)
				{
					int jcol = LocExt_colind[j];

					if (LocExt_s[j] && Ext_CF_map[jcol] == GREATER)
					{
						if (Ext_weights[jcol] > weights[irow])
						{
							CF_map[irow] = LESS;
							break;
						}
					}

				}
			}
		}

		for (int r = 0; r < num_neighbors; ++r)
			if (recv_ptr[r + 1] > recv_ptr[r] && neighbor_unassigned[r] != -1)
				MPI_Irecv(Ext_CF_map + recv_ptr[r], recv_ptr[r + 1] - recv_ptr[r], MPI_INT, neighbor_ranks[r], MPI_TAG, comm, recv_requests + r);

		for (int r = 0; r < num_neighbors; ++r)
		{
			if (send_ptr[r + 1] > send_ptr[r] && neighbor_unassigned[r] != -1)
			{
				for (int i = send_ptr[r]; i < send_ptr[r + 1]; ++i)
					send_buffer[i] = CF_map[send_list[i]];
				MPI_Isend(send_buffer + send_ptr[r], send_ptr[r + 1] - send_ptr[r], MPI_INT, neighbor_ranks[r], MPI_TAG, comm, send_requests + r);
			}
		}

		for (int r = 0; r < num_neighbors; ++r)
			if (recv_ptr[r + 1] > recv_ptr[r] && neighbor_unassigned[r] != -1)
				MPI_Wait(recv_requests + r, &status);

		for (int r = 0; r < num_neighbors; ++r)
			if (send_ptr[r + 1] > send_ptr[r] && neighbor_unassigned[r] != -1)
				MPI_Wait(send_requests + r, &status);

#ifdef USE_OPENMP
#pragma omp parallel for reduction(+:num_coarse) schedule(guided) 
#endif
		for (int i = 0; i < n; ++i)
		{
			if (CF_map[i] == GREATER)
			{
				CF_map[i] = COARSE;
				++num_coarse;
			}
		}

#ifdef USE_OPENMP
#pragma omp parallel for reduction(+:num_fine) schedule(guided) 
#endif
		for (int i = 0; i < n; ++i)
		{
			if (CF_map[i] == LESS)
			{
				for (int j = LocLoc_rowptr[i]; j < LocLoc_rowptr[i + 1]; ++j)
				{
					if (LocLoc_s[j] && CF_map[LocLoc_colind[j]] == COARSE)
					{
						CF_map[i] = FINE;
						++num_fine;
						break;
					}
				}
			}
		}

#ifdef USE_OPENMP
#pragma omp parallel for reduction(+:num_fine) schedule(guided)
#endif
		for (int i = 0; i < num_interface_rows; ++i)
		{
			int irow = interface_rows[i];

			if (CF_map[irow] == LESS)
			{
				for (int j = LocExt_rowptr[irow]; j < LocExt_rowptr[irow + 1]; ++j)
				{
					if (LocExt_s[j] && Ext_CF_map[LocExt_colind[j]] == GREATER)
					{
						CF_map[irow] = FINE;
						++num_fine;
						break;
					}
				}
			}
		}

		for (int r = 0; r < num_neighbors; ++r)
			if (recv_ptr[r + 1] > recv_ptr[r] && neighbor_unassigned[r] != -1)
				MPI_Irecv(Ext_CF_map + recv_ptr[r], recv_ptr[r + 1] - recv_ptr[r], MPI_INT, neighbor_ranks[r], MPI_TAG, comm, recv_requests + r);

		for (int r = 0; r < num_neighbors; ++r)
		{
			if (send_ptr[r + 1] > send_ptr[r] && neighbor_unassigned[r] != -1)
			{
				for (int i = send_ptr[r]; i < send_ptr[r + 1]; ++i)
					send_buffer[i] = CF_map[send_list[i]];
				MPI_Isend(send_buffer + send_ptr[r], send_ptr[r + 1] - send_ptr[r], MPI_INT, neighbor_ranks[r], MPI_TAG, comm, send_requests + r);
			}
		}

		for (int r = 0; r < num_neighbors; ++r)
			if (recv_ptr[r + 1] > recv_ptr[r] && neighbor_unassigned[r] != -1)
				MPI_Wait(recv_requests + r, &status);

		for (int r = 0; r < num_neighbors; ++r)
			if (send_ptr[r + 1] > send_ptr[r] && neighbor_unassigned[r] != -1)
				MPI_Wait(send_requests + r, &status);
	}

	delete[] recv_requests;
	delete[] send_requests;
	delete[] send_buffer;
	delete[] weights;
	delete[] interface_rows;
	delete[] Ext_weights;
	delete[] Ext_CF_map;
	delete[] neighbor_unassigned;

	return num_coarse;
}

static int ParCSR_HMIS_Coarsening(const ParCSR_Matrix& S2, int* CF_map)
{
	MPI_Comm comm = S2.comm;

	int rank;
	MPI_Comm_rank(comm, &rank);

	int n = S2.LocLoc.size[0];
	int Ext_size = S2.LocExt.size[1];

	int* LocLoc_rowptr = S2.LocLoc.rowptr;
	int* LocLoc_colind = S2.LocLoc.colind;
	int* LocExt_rowptr = S2.LocExt.rowptr;
	int* LocExt_colind = S2.LocExt.colind;

	int num_interface_rows = 0;
	int* interface_rows = new int[n];
	for (int i = 0; i < n; ++i)
		if (LocExt_rowptr[i + 1] > LocExt_rowptr[i])
			interface_rows[num_interface_rows++] = i;

	int num_neighbors = S2.num_neighbors;
	int* neighbor_ranks = S2.neighbor_ranks;
	int* recv_ptr = S2.Ext_ptr;
	int* send_ptr = S2.send_ptr;
	int* send_list = S2.send_list;
	int* send_buffer = new int[send_ptr[num_neighbors]];

	int* weights = new int[n];
	int* Ext_weights = new int[Ext_size];

	int* S2tp = new int[n + 1];
	for (int i = 0; i <= n; ++i)
		S2tp[i] = 0;

	for (int j = LocLoc_rowptr[0]; j < LocLoc_rowptr[n]; ++j)
		++S2tp[LocLoc_colind[j]];

	for (int i = 0; i < n; ++i)
		S2tp[i + 1] += S2tp[i];

	int* S2ti = new int[S2tp[n]];

	for (int i = n - 1; i >= 0; --i)
		for (int j = LocLoc_rowptr[i + 1] - 1; j >= LocLoc_rowptr[i]; --j)
			S2ti[--S2tp[LocLoc_colind[j]]] = i;

	for (int i = 0; i < n; ++i)
		weights[i] = S2tp[i + 1] - S2tp[i];

	int num_coarse = 0;
	int num_fine = 0;

	for (int i = 0; i < n; ++i)
		CF_map[i] = UNASSIGNED;

	while (num_coarse + num_fine < n)
	{
		for (int i = 0; i < n; ++i)
		{
			int r = i;

			if (CF_map[r] == UNASSIGNED)
			{
			L0:
				for (int j = LocLoc_rowptr[r]; j < LocLoc_rowptr[r + 1]; ++j)
				{
					int jcol = LocLoc_colind[j];

					if (CF_map[jcol] == UNASSIGNED)
					{
						if (weights[r] < weights[jcol])
						{
							r = jcol;
							goto L0;
						}
					}
				}

				for (int j = S2tp[r]; j < S2tp[r + 1]; ++j)
				{
					int jcol = S2ti[j];

					if (CF_map[jcol] == UNASSIGNED)
					{
						if (weights[r] < weights[jcol])
						{
							r = jcol;
							goto L0;
						}
					}
				}

				CF_map[r] = COARSE;
				++num_coarse;

				for (int j = S2tp[r]; j < S2tp[r + 1]; ++j)
				{
					int jcol = S2ti[j];
					if (CF_map[jcol] == UNASSIGNED)
					{
						CF_map[jcol] = FINE;
						++num_fine;

						for (int k = LocLoc_rowptr[jcol]; k < LocLoc_rowptr[jcol + 1]; ++k)
						{
							int kcol = LocLoc_colind[k];

							if (CF_map[kcol] == UNASSIGNED)
								++weights[kcol];
						}
					}
				}

				for (int j = LocLoc_rowptr[r]; j < LocLoc_rowptr[r + 1]; ++j)
				{
					int jcol = LocLoc_colind[j];

					if (CF_map[jcol] == UNASSIGNED)
					{
						if (--weights[jcol] == 0)
						{
							CF_map[jcol] = FINE;
							++num_fine;

							for (int k = LocLoc_rowptr[jcol]; k < LocLoc_rowptr[jcol + 1]; ++k)
							{
								int kcol = LocLoc_colind[k];

								if (CF_map[kcol] == UNASSIGNED)
									++weights[kcol];
							}
						}
					}
				}
			}
		}
	}

	int lambda = RANDOM_MAX + 1;
	static unsigned int x = rank + 2671;

	for (int i = 0; i < num_interface_rows; ++i)
	{
		int irow = interface_rows[i];
		weights[irow] = (S2tp[irow + 1] - S2tp[irow]) * lambda + Random(x);
	}

	for (int i = 0; i < Ext_size; ++i)
		Ext_weights[i] = 0;

	for (int j = LocExt_rowptr[0]; j < LocExt_rowptr[n]; ++j)
		Ext_weights[LocExt_colind[j]] += lambda;

	MPI_Request* recv_requests = new MPI_Request[num_neighbors];
	MPI_Request* send_requests = new MPI_Request[num_neighbors];
	MPI_Status status;

	for (int r = 0; r < num_neighbors; ++r)
		if (send_ptr[r + 1] > send_ptr[r])
			MPI_Irecv(send_buffer + send_ptr[r], send_ptr[r + 1] - send_ptr[r], MPI_INT, neighbor_ranks[r], MPI_TAG, comm, recv_requests + r);

	for (int r = 0; r < num_neighbors; ++r)
		if (recv_ptr[r + 1] > recv_ptr[r])
			MPI_Isend(Ext_weights + recv_ptr[r], recv_ptr[r + 1] - recv_ptr[r], MPI_INT, neighbor_ranks[r], MPI_TAG, comm, send_requests + r);

	for (int r = 0; r < num_neighbors; ++r)
	{
		if (send_ptr[r + 1] > send_ptr[r])
		{
			MPI_Wait(recv_requests + r, &status);
			for (int i = send_ptr[r]; i < send_ptr[r + 1]; ++i)
				weights[send_list[i]] += send_buffer[i];
		}
	}

	for (int r = 0; r < num_neighbors; ++r)
		if (recv_ptr[r + 1] > recv_ptr[r])
			MPI_Wait(send_requests + r, &status);

	for (int r = 0; r < num_neighbors; ++r)
		if (recv_ptr[r + 1] > recv_ptr[r])
			MPI_Irecv(Ext_weights + recv_ptr[r], recv_ptr[r + 1] - recv_ptr[r], MPI_INT, neighbor_ranks[r], MPI_TAG, comm, recv_requests + r);

	for (int r = 0; r < num_neighbors; ++r)
	{
		if (send_ptr[r + 1] > send_ptr[r])
		{
			for (int i = send_ptr[r]; i < send_ptr[r + 1]; ++i)
				send_buffer[i] = weights[send_list[i]];
			MPI_Isend(send_buffer + send_ptr[r], send_ptr[r + 1] - send_ptr[r], MPI_INT, neighbor_ranks[r], MPI_TAG, comm, send_requests + r);
		}
	}

	for (int r = 0; r < num_neighbors; ++r)
		if (recv_ptr[r + 1] > recv_ptr[r])
			MPI_Wait(recv_requests + r, &status);

	for (int r = 0; r < num_neighbors; ++r)
		if (send_ptr[r + 1] > send_ptr[r])
			MPI_Wait(send_requests + r, &status);

	delete[] S2tp;
	delete[] S2ti;

	for (int i = 0; i < num_interface_rows; ++i)
	{
		int irow = interface_rows[i];
		--(CF_map[irow] == COARSE ? num_coarse : num_fine);
		CF_map[irow] = UNASSIGNED;
	}

	for (int i = 0; i < num_interface_rows; ++i)
	{
		int irow = interface_rows[i];

		for (int j = LocLoc_rowptr[irow]; j < LocLoc_rowptr[irow + 1]; ++j)
		{
			int jcol = LocLoc_colind[j];

			if (CF_map[jcol] == COARSE)
			{
				CF_map[irow] = FINE;
				++num_fine;
				break;
			}
		}
	}

	int* Ext_CF_map = new int[Ext_size];

	for (int r = 0; r < num_neighbors; ++r)
		if (recv_ptr[r + 1] > recv_ptr[r])
			MPI_Irecv(Ext_CF_map + recv_ptr[r], recv_ptr[r + 1] - recv_ptr[r], MPI_INT, neighbor_ranks[r], MPI_TAG, comm, recv_requests + r);

	for (int r = 0; r < num_neighbors; ++r)
	{
		if (send_ptr[r + 1] > send_ptr[r])
		{
			for (int i = send_ptr[r]; i < send_ptr[r + 1]; ++i)
				send_buffer[i] = CF_map[send_list[i]];
			MPI_Isend(send_buffer + send_ptr[r], send_ptr[r + 1] - send_ptr[r], MPI_INT, neighbor_ranks[r], MPI_TAG, comm, send_requests + r);
		}
	}

	for (int r = 0; r < num_neighbors; ++r)
		if (recv_ptr[r + 1] > recv_ptr[r])
			MPI_Wait(recv_requests + r, &status);

	for (int r = 0; r < num_neighbors; ++r)
		if (send_ptr[r + 1] > send_ptr[r])
			MPI_Wait(send_requests + r, &status);

	int* neighbor_unassigned = new int[num_neighbors];
	for (int r = 0; r < num_neighbors; ++r)
		neighbor_unassigned[r] = 1;

	int LESS = 0;
	int GREATER = 2;

	while (1)
	{
		int num_unassigned = n - num_coarse - num_fine;

		for (int r = 0; r < num_neighbors; ++r)
			if (neighbor_unassigned[r] != -1)
				MPI_Irecv(neighbor_unassigned + r, 1, MPI_INT, neighbor_ranks[r], MPI_TAG, comm, recv_requests + r);

		for (int r = 0; r < num_neighbors; ++r)
			if (neighbor_unassigned[r] != -1)
				MPI_Send(&num_unassigned, 1, MPI_INT, neighbor_ranks[r], MPI_TAG, comm);

		for (int r = 0; r < num_neighbors; ++r)
		{
			if (neighbor_unassigned[r] != -1)
			{
				MPI_Wait(recv_requests + r, &status);
				if (neighbor_unassigned[r] == 0)
					neighbor_unassigned[r] = -1;
			}
		}

		if (num_unassigned == 0) break;

		Swap(LESS, GREATER);

		for (int i = 0; i < num_interface_rows; ++i)
		{
			int irow = interface_rows[i];

			if (CF_map[irow] == GREATER)
			{
				for (int j = LocLoc_rowptr[irow]; j < LocLoc_rowptr[irow + 1]; ++j)
				{
					int jcol = LocLoc_colind[j];

					if (CF_map[jcol] == GREATER)
					{
						if (weights[jcol] > weights[irow])
						{
							CF_map[irow] = LESS;
							break;
						}
						else
							CF_map[jcol] = LESS;
					}
				}
			}
		}

		for (int i = 0; i < num_interface_rows; ++i)
		{
			int irow = interface_rows[i];

			if (CF_map[irow] == GREATER)
			{
				for (int j = LocExt_rowptr[irow]; j < LocExt_rowptr[irow + 1]; ++j)
				{
					int jcol = LocExt_colind[j];

					if (Ext_CF_map[jcol] == GREATER)
					{
						if (Ext_weights[jcol] > weights[irow])
						{
							CF_map[irow] = LESS;
							break;
						}
					}

				}
			}
		}

		for (int r = 0; r < num_neighbors; ++r)
			if (recv_ptr[r + 1] > recv_ptr[r] && neighbor_unassigned[r] != -1)
				MPI_Irecv(Ext_CF_map + recv_ptr[r], recv_ptr[r + 1] - recv_ptr[r], MPI_INT, neighbor_ranks[r], MPI_TAG, comm, recv_requests + r);

		for (int r = 0; r < num_neighbors; ++r)
		{
			if (send_ptr[r + 1] > send_ptr[r] && neighbor_unassigned[r] != -1)
			{
				for (int i = send_ptr[r]; i < send_ptr[r + 1]; ++i)
					send_buffer[i] = CF_map[send_list[i]];
				MPI_Isend(send_buffer + send_ptr[r], send_ptr[r + 1] - send_ptr[r], MPI_INT, neighbor_ranks[r], MPI_TAG, comm, send_requests + r);
			}
		}

		for (int r = 0; r < num_neighbors; ++r)
			if (recv_ptr[r + 1] > recv_ptr[r] && neighbor_unassigned[r] != -1)
				MPI_Wait(recv_requests + r, &status);

		for (int r = 0; r < num_neighbors; ++r)
			if (send_ptr[r + 1] > send_ptr[r] && neighbor_unassigned[r] != -1)
				MPI_Wait(send_requests + r, &status);

		for (int i = 0; i < num_interface_rows; ++i)
		{
			int irow = interface_rows[i];

			if (CF_map[irow] == GREATER)
			{
				CF_map[irow] = COARSE;
				++num_coarse;
			}
		}

		for (int i = 0; i < num_interface_rows; ++i)
		{
			int irow = interface_rows[i];

			if (CF_map[irow] == LESS)
			{
				for (int j = LocLoc_rowptr[irow]; j < LocLoc_rowptr[irow + 1]; ++j)
				{
					if (CF_map[LocLoc_colind[j]] == COARSE)
					{
						CF_map[irow] = FINE;
						++num_fine;
						break;
					}
				}
			}
		}

		for (int i = 0; i < num_interface_rows; ++i)
		{
			int irow = interface_rows[i];

			if (CF_map[irow] == LESS)
			{
				for (int j = LocExt_rowptr[irow]; j < LocExt_rowptr[irow + 1]; ++j)
				{
					if (Ext_CF_map[LocExt_colind[j]] == GREATER)
					{
						CF_map[irow] = FINE;
						++num_fine;
						break;
					}
				}
			}
		}

		for (int r = 0; r < num_neighbors; ++r)
			if (recv_ptr[r + 1] > recv_ptr[r] && neighbor_unassigned[r] != -1)
				MPI_Irecv(Ext_CF_map + recv_ptr[r], recv_ptr[r + 1] - recv_ptr[r], MPI_INT, neighbor_ranks[r], MPI_TAG, comm, recv_requests + r);

		for (int r = 0; r < num_neighbors; ++r)
		{
			if (send_ptr[r + 1] > send_ptr[r] && neighbor_unassigned[r] != -1)
			{
				for (int i = send_ptr[r]; i < send_ptr[r + 1]; ++i)
					send_buffer[i] = CF_map[send_list[i]];
				MPI_Isend(send_buffer + send_ptr[r], send_ptr[r + 1] - send_ptr[r], MPI_INT, neighbor_ranks[r], MPI_TAG, comm, send_requests + r);
			}
		}

		for (int r = 0; r < num_neighbors; ++r)
			if (recv_ptr[r + 1] > recv_ptr[r] && neighbor_unassigned[r] != -1)
				MPI_Wait(recv_requests + r, &status);

		for (int r = 0; r < num_neighbors; ++r)
			if (send_ptr[r + 1] > send_ptr[r] && neighbor_unassigned[r] != -1)
				MPI_Wait(send_requests + r, &status);
	}

	delete[] recv_requests;
	delete[] send_requests;
	delete[] send_buffer;
	delete[] weights;
	delete[] interface_rows;
	delete[] Ext_weights;
	delete[] Ext_CF_map;
	delete[] neighbor_unassigned;

	return num_coarse;
}

static int ParCSR_PMIS_Coarsening(const ParCSR_Matrix& S2, int* CF_map)
{
	MPI_Comm comm = S2.comm;

	int rank;
	MPI_Comm_rank(comm, &rank);

	int n = S2.LocLoc.size[0];
	int Ext_size = S2.LocExt.size[1];

	int* LocLoc_rowptr = S2.LocLoc.rowptr;
	int* LocLoc_colind = S2.LocLoc.colind;
	int* LocExt_rowptr = S2.LocExt.rowptr;
	int* LocExt_colind = S2.LocExt.colind;

	int num_interface_rows = 0;
	int* interface_rows = new int[n];
	for (int i = 0; i < n; ++i)
		if (LocExt_rowptr[i + 1] > LocExt_rowptr[i])
			interface_rows[num_interface_rows++] = i;

	int num_neighbors = S2.num_neighbors;
	int* neighbor_ranks = S2.neighbor_ranks;
	int* recv_ptr = S2.Ext_ptr;
	int* send_ptr = S2.send_ptr;
	int* send_list = S2.send_list;
	int* send_buffer = new int[send_ptr[num_neighbors]];

	int* weights = new int[n];
	int* Ext_weights = new int[Ext_size];

	MPI_Request* recv_requests = new MPI_Request[num_neighbors];
	MPI_Request* send_requests = new MPI_Request[num_neighbors];
	MPI_Status status;

	int lambda = RANDOM_MAX + 1;
	static unsigned int x = rank + 2671;

	for (int i = 0; i < n; ++i)
		weights[i] = Random(x);

	for (int i = 0; i < Ext_size; ++i)
		Ext_weights[i] = 0;

	for (int j = LocLoc_rowptr[0]; j < LocLoc_rowptr[n]; ++j)
		weights[LocLoc_colind[j]] += lambda;

	for (int j = LocExt_rowptr[0]; j < LocExt_rowptr[n]; ++j)
		Ext_weights[LocExt_colind[j]] += lambda;

	for (int r = 0; r < num_neighbors; ++r)
		if (send_ptr[r + 1] > send_ptr[r])
			MPI_Irecv(send_buffer + send_ptr[r], send_ptr[r + 1] - send_ptr[r], MPI_INT, neighbor_ranks[r], MPI_TAG, comm, recv_requests + r);

	for (int r = 0; r < num_neighbors; ++r)
		if (recv_ptr[r + 1] > recv_ptr[r])
			MPI_Isend(Ext_weights + recv_ptr[r], recv_ptr[r + 1] - recv_ptr[r], MPI_INT, neighbor_ranks[r], MPI_TAG, comm, send_requests + r);

	for (int r = 0; r < num_neighbors; ++r)
	{
		if (send_ptr[r + 1] > send_ptr[r])
		{
			MPI_Wait(recv_requests + r, &status);
			for (int i = send_ptr[r]; i < send_ptr[r + 1]; ++i)
				weights[send_list[i]] += send_buffer[i];
		}
	}

	for (int r = 0; r < num_neighbors; ++r)
		if (recv_ptr[r + 1] > recv_ptr[r])
			MPI_Wait(send_requests + r, &status);

	for (int r = 0; r < num_neighbors; ++r)
		if (recv_ptr[r + 1] > recv_ptr[r])
			MPI_Irecv(Ext_weights + recv_ptr[r], recv_ptr[r + 1] - recv_ptr[r], MPI_INT, neighbor_ranks[r], MPI_TAG, comm, recv_requests + r);

	for (int r = 0; r < num_neighbors; ++r)
	{
		if (send_ptr[r + 1] > send_ptr[r])
		{
			for (int i = send_ptr[r]; i < send_ptr[r + 1]; ++i)
				send_buffer[i] = weights[send_list[i]];
			MPI_Isend(send_buffer + send_ptr[r], send_ptr[r + 1] - send_ptr[r], MPI_INT, neighbor_ranks[r], MPI_TAG, comm, send_requests + r);
		}
	}

	for (int r = 0; r < num_neighbors; ++r)
		if (recv_ptr[r + 1] > recv_ptr[r])
			MPI_Wait(recv_requests + r, &status);

	for (int r = 0; r < num_neighbors; ++r)
		if (send_ptr[r + 1] > send_ptr[r])
			MPI_Wait(send_requests + r, &status);

	int num_coarse = 0;
	int num_fine = 0;

	for (int i = 0; i < n; ++i)
		CF_map[i] = UNASSIGNED;

	int* Ext_CF_map = new int[Ext_size];
	for (int i = 0; i < Ext_size; ++i)
		Ext_CF_map[i] = UNASSIGNED;

	int* neighbor_unassigned = new int[num_neighbors];
	for (int r = 0; r < num_neighbors; ++r)
		neighbor_unassigned[r] = 1;

	int LESS = 0;
	int GREATER = 2;

	while (1)
	{
		int num_unassigned = n - num_coarse - num_fine;

		for (int r = 0; r < num_neighbors; ++r)
			if (neighbor_unassigned[r] != -1)
				MPI_Irecv(neighbor_unassigned + r, 1, MPI_INT, neighbor_ranks[r], MPI_TAG, comm, recv_requests + r);

		for (int r = 0; r < num_neighbors; ++r)
			if (neighbor_unassigned[r] != -1)
				MPI_Send(&num_unassigned, 1, MPI_INT, neighbor_ranks[r], MPI_TAG, comm);

		for (int r = 0; r < num_neighbors; ++r)
		{
			if (neighbor_unassigned[r] != -1)
			{
				MPI_Wait(recv_requests + r, &status);
				if (neighbor_unassigned[r] == 0)
					neighbor_unassigned[r] = -1;
			}
		}

		if (num_unassigned == 0) break;

		Swap(LESS, GREATER);

#ifdef USE_OPENMP
#pragma omp parallel for schedule(guided)
#endif
		for (int i = 0; i < n; ++i)
		{
			if (CF_map[i] == GREATER)
			{
				for (int j = LocLoc_rowptr[i]; j < LocLoc_rowptr[i + 1]; ++j)
				{
					int jcol = LocLoc_colind[j];

					if (CF_map[jcol] == GREATER)
					{
						if (weights[jcol] > weights[i])
						{
							CF_map[i] = LESS;
							break;
						}
						else
							CF_map[jcol] = LESS;
					}
				}
			}
		}

#ifdef USE_OPENMP
#pragma omp parallel for schedule(guided)
#endif
		for (int i = 0; i < num_interface_rows; ++i)
		{
			int irow = interface_rows[i];

			if (CF_map[irow] == GREATER)
			{
				for (int j = LocExt_rowptr[irow]; j < LocExt_rowptr[irow + 1]; ++j)
				{
					int jcol = LocExt_colind[j];

					if (Ext_CF_map[jcol] == GREATER)
					{
						if (Ext_weights[jcol] > weights[irow])
						{
							CF_map[irow] = LESS;
							break;
						}
					}
				}
			}
		}

		for (int r = 0; r < num_neighbors; ++r)
			if (recv_ptr[r + 1] > recv_ptr[r] && neighbor_unassigned[r] != -1)
				MPI_Irecv(Ext_CF_map + recv_ptr[r], recv_ptr[r + 1] - recv_ptr[r], MPI_INT, neighbor_ranks[r], MPI_TAG, comm, recv_requests + r);

		for (int r = 0; r < num_neighbors; ++r)
		{
			if (send_ptr[r + 1] > send_ptr[r] && neighbor_unassigned[r] != -1)
			{
				for (int i = send_ptr[r]; i < send_ptr[r + 1]; ++i)
					send_buffer[i] = CF_map[send_list[i]];
				MPI_Isend(send_buffer + send_ptr[r], send_ptr[r + 1] - send_ptr[r], MPI_INT, neighbor_ranks[r], MPI_TAG, comm, send_requests + r);
			}
		}

		for (int r = 0; r < num_neighbors; ++r)
			if (recv_ptr[r + 1] > recv_ptr[r] && neighbor_unassigned[r] != -1)
				MPI_Wait(recv_requests + r, &status);

		for (int r = 0; r < num_neighbors; ++r)
			if (send_ptr[r + 1] > send_ptr[r] && neighbor_unassigned[r] != -1)
				MPI_Wait(send_requests + r, &status);

#ifdef USE_OPENMP
#pragma omp parallel for reduction(+:num_coarse) schedule(guided) 
#endif
		for (int i = 0; i < n; ++i)
		{
			if (CF_map[i] == GREATER)
			{
				CF_map[i] = COARSE;
				++num_coarse;
			}
		}

#ifdef USE_OPENMP
#pragma omp parallel for reduction(+:num_fine) schedule(guided) 
#endif
		for (int i = 0; i < n; ++i)
		{
			if (CF_map[i] == LESS)
			{
				for (int j = LocLoc_rowptr[i]; j < LocLoc_rowptr[i + 1]; ++j)
				{
					if (CF_map[LocLoc_colind[j]] == COARSE)
					{
						CF_map[i] = FINE;
						++num_fine;
						break;
					}
				}
			}
		}

#ifdef USE_OPENMP
#pragma omp parallel for reduction(+:num_fine) schedule(guided)
#endif
		for (int i = 0; i < num_interface_rows; ++i)
		{
			int irow = interface_rows[i];

			if (CF_map[irow] == LESS)
			{
				for (int j = LocExt_rowptr[irow]; j < LocExt_rowptr[irow + 1]; ++j)
				{
					if (Ext_CF_map[LocExt_colind[j]] == GREATER)
					{
						CF_map[irow] = FINE;
						++num_fine;
						break;
					}
				}
			}
		}

		for (int r = 0; r < num_neighbors; ++r)
			if (recv_ptr[r + 1] > recv_ptr[r] && neighbor_unassigned[r] != -1)
				MPI_Irecv(Ext_CF_map + recv_ptr[r], recv_ptr[r + 1] - recv_ptr[r],
					MPI_INT, neighbor_ranks[r], MPI_TAG, comm, recv_requests + r);

		for (int r = 0; r < num_neighbors; ++r)
		{
			if (send_ptr[r + 1] > send_ptr[r] && neighbor_unassigned[r] != -1)
			{
				for (int i = send_ptr[r]; i < send_ptr[r + 1]; ++i)
					send_buffer[i] = CF_map[send_list[i]];
				MPI_Isend(send_buffer + send_ptr[r], send_ptr[r + 1] - send_ptr[r],
					MPI_INT, neighbor_ranks[r], MPI_TAG, comm, send_requests + r);
			}
		}

		for (int r = 0; r < num_neighbors; ++r)
			if (recv_ptr[r + 1] > recv_ptr[r] && neighbor_unassigned[r] != -1)
				MPI_Wait(recv_requests + r, &status);

		for (int r = 0; r < num_neighbors; ++r)
			if (send_ptr[r + 1] > send_ptr[r] && neighbor_unassigned[r] != -1)
				MPI_Wait(send_requests + r, &status);
	}

	delete[] recv_requests;
	delete[] send_requests;
	delete[] send_buffer;
	delete[] weights;
	delete[] interface_rows;
	delete[] Ext_weights;
	delete[] Ext_CF_map;
	delete[] neighbor_unassigned;

	return num_coarse;
}

static void ParCSR_Strength2(const ParCSR_Matrix& A, const bool* LocLoc_s, const bool* LocExt_s, int* CF_map, ParCSR_Matrix& S2)
{
	MPI_Comm comm = A.comm;

	int rank;
	MPI_Comm_rank(comm, &rank);

	int n = A.LocLoc.size[0];
	int Ext_size = A.LocExt.size[1];

	int* A_LocLoc_rowptr = A.LocLoc.rowptr;
	int* A_LocLoc_colind = A.LocLoc.colind;
	int* A_LocExt_rowptr = A.LocExt.rowptr;
	int* A_LocExt_colind = A.LocExt.colind;

	int num_neighbors = A.num_neighbors;
	int* neighbor_ranks = A.neighbor_ranks;
	int* recv_ptr = A.Ext_ptr;
	int* send_ptr = A.send_ptr;
	int* send_list = A.send_list;
	int* send_cf_m = new int[send_ptr[num_neighbors]];

	int num_coarse = 0;
	for (int i = 0; i < n; ++i)
		if (CF_map[i] != FINE)
			CF_map[i] = num_coarse++;

	int* Ext_CF_map = new int[Ext_size];

	MPI_Request* recv_requests = new MPI_Request[num_neighbors];
	MPI_Request* send_requests = new MPI_Request[num_neighbors];
	MPI_Status status;

	for (int r = 0; r < num_neighbors; ++r)
		if (recv_ptr[r + 1] > recv_ptr[r])
			MPI_Irecv(Ext_CF_map + recv_ptr[r], recv_ptr[r + 1] - recv_ptr[r], MPI_INT, neighbor_ranks[r], MPI_TAG, comm, recv_requests + r);

	for (int r = 0; r < num_neighbors; ++r)
	{
		if (send_ptr[r + 1] > send_ptr[r])
		{
			for (int i = send_ptr[r]; i < send_ptr[r + 1]; ++i)
				send_cf_m[i] = CF_map[send_list[i]];
			MPI_Isend(send_cf_m + send_ptr[r], send_ptr[r + 1] - send_ptr[r], MPI_INT, neighbor_ranks[r], MPI_TAG, comm, send_requests + r);
		}
	}

	for (int r = 0; r < num_neighbors; ++r)
		if (recv_ptr[r + 1] > recv_ptr[r])
			MPI_Wait(recv_requests + r, &status);

	for (int r = 0; r < num_neighbors; ++r)
		if (send_ptr[r + 1] > send_ptr[r])
			MPI_Wait(send_requests + r, &status);

	delete[] send_cf_m;

	int* Sc_Ext_map = new int[Ext_size];
	int* Sc_Ext_ptr = new int[num_neighbors + 1];

	Sc_Ext_ptr[0] = 0;
	for (int r = 0, k = 0; r < num_neighbors; Sc_Ext_ptr[++r] = k)
	{
		for (int i = recv_ptr[r]; i < recv_ptr[r + 1]; ++i)
		{
			if (Ext_CF_map[i] != FINE)
			{
				Sc_Ext_map[k] = Ext_CF_map[i];
				Ext_CF_map[i] = k++;
			}
		}
	}

	int* Sc_LocLoc_rowptr = new int[n + 1];
	int* Sc_LocExt_rowptr = new int[n + 1];

	Sc_LocLoc_rowptr[0] = 0;
	Sc_LocExt_rowptr[0] = 0;

#ifdef USE_OPENMP
#pragma omp parallel for schedule(guided)
#endif
	for (int i = 0; i < n; ++i)
	{
		if (CF_map[i] == FINE)
		{
			int cnt_loc = 0;
			int cnt_ext = 0;

			for (int j = A_LocLoc_rowptr[i]; j < A_LocLoc_rowptr[i + 1]; ++j)
				if (LocLoc_s[j] && CF_map[A_LocLoc_colind[j]] != FINE) ++cnt_loc;
			for (int j = A_LocExt_rowptr[i]; j < A_LocExt_rowptr[i + 1]; ++j)
				if (LocExt_s[j] && Ext_CF_map[A_LocExt_colind[j]] != FINE) ++cnt_ext;

			Sc_LocLoc_rowptr[i + 1] = cnt_loc;
			Sc_LocExt_rowptr[i + 1] = cnt_ext;
		}
		else
		{
			Sc_LocLoc_rowptr[i + 1] = 0;
			Sc_LocExt_rowptr[i + 1] = 0;
		}
	}

	for (int i = 0; i < n; ++i)
		Sc_LocLoc_rowptr[i + 1] += Sc_LocLoc_rowptr[i];
	for (int i = 0; i < n; ++i)
		Sc_LocExt_rowptr[i + 1] += Sc_LocExt_rowptr[i];

	int* Sc_LocLoc_colind = new int[Sc_LocLoc_rowptr[n]];
	int* Sc_LocExt_colind = new int[Sc_LocExt_rowptr[n]];

#ifdef USE_OPENMP
#pragma omp parallel for schedule(guided)
#endif
	for (int i = 0; i < n; ++i)
	{
		if (CF_map[i] == FINE)
		{
			for (int j = A_LocLoc_rowptr[i], k = Sc_LocLoc_rowptr[i]; j < A_LocLoc_rowptr[i + 1]; ++j)
				if (LocLoc_s[j] && CF_map[A_LocLoc_colind[j]] != FINE)
					Sc_LocLoc_colind[k++] = CF_map[A_LocLoc_colind[j]];

			for (int j = A_LocExt_rowptr[i], k = Sc_LocExt_rowptr[i]; j < A_LocExt_rowptr[i + 1]; ++j)
				if (LocExt_s[j] && Ext_CF_map[A_LocExt_colind[j]] != FINE)
					Sc_LocExt_colind[k++] = Ext_CF_map[A_LocExt_colind[j]];
		}
	}

	int* Sc_Ext_from = new int[Sc_Ext_ptr[num_neighbors]];
	for (int r = 0; r < num_neighbors; ++r)
		for (int i = Sc_Ext_ptr[r]; i < Sc_Ext_ptr[r + 1]; ++i)
			Sc_Ext_from[i] = r;

	std::map<int, int> map_neighbor;
	for (int r = 0; r < num_neighbors; ++r)
		map_neighbor.insert(std::pair<int, int>(neighbor_ranks[r], r));

	std::set<int>* set_neighbors = new std::set<int>[num_neighbors];

	for (int r = 0; r < num_neighbors; ++r)
	{
		for (int i = send_ptr[r]; i < send_ptr[r + 1]; ++i)
		{
			int irow = send_list[i];
			for (int j = Sc_LocExt_rowptr[irow]; j < Sc_LocExt_rowptr[irow + 1]; ++j)
			{
				int jfrom = Sc_Ext_from[Sc_LocExt_colind[j]];
				if (jfrom != r)
				{
					set_neighbors[r].insert(neighbor_ranks[jfrom]);
					set_neighbors[jfrom].insert(neighbor_ranks[r]);
				}
			}
		}
	}

	int** recv_s2_nb = new int* [num_neighbors];
	int** send_s2_nb = new int* [num_neighbors];

	for (int r = 0; r < num_neighbors; ++r)
	{
		send_s2_nb[r] = new int[set_neighbors[r].size()];
		int i = 0;
		for (std::set<int>::iterator iter = set_neighbors[r].begin(); iter != set_neighbors[r].end(); ++iter)
			send_s2_nb[r][i++] = *iter;
		MPI_Isend(send_s2_nb[r], set_neighbors[r].size(), MPI_INT, neighbor_ranks[r], MPI_TAG, comm, send_requests + r);
	}

	for (int r = 0; r < num_neighbors; ++r)
	{
		int cnt;
		MPI_Probe(neighbor_ranks[r], MPI_TAG, comm, &status);
		MPI_Get_count(&status, MPI_INT, &cnt);
		recv_s2_nb[r] = new int[cnt];
		MPI_Irecv(recv_s2_nb[r], cnt, MPI_INT, neighbor_ranks[r], MPI_TAG, comm, recv_requests + r);
	}

	for (int r = 0; r < num_neighbors; ++r)
	{
		int cnt;
		MPI_Wait(recv_requests + r, &status);
		MPI_Get_count(&status, MPI_INT, &cnt);
		for (int i = 0; i < cnt; ++i)
			map_neighbor.insert(std::pair<int, int>(recv_s2_nb[r][i], map_neighbor.size()));
		delete[] recv_s2_nb[r];
	}

	for (int r = 0; r < num_neighbors; ++r)
	{
		MPI_Wait(send_requests + r, &status);
		delete[] send_s2_nb[r];
	}

	delete[] set_neighbors;
	delete[] recv_s2_nb;
	delete[] send_s2_nb;

	int S2_num_neighbors = map_neighbor.size();
	int* S2_neighbors = new int[S2_num_neighbors];
	for (std::map<int, int>::iterator iter = map_neighbor.begin(); iter != map_neighbor.end(); ++iter)
		S2_neighbors[iter->second] = iter->first;

	int* send_bufsize = new int[num_neighbors];
	int* recv_bufsize = new int[num_neighbors];

	for (int r = 0; r < num_neighbors; ++r)
		if (recv_ptr[r + 1] > recv_ptr[r])
			MPI_Irecv(recv_bufsize + r, 1, MPI_INT, neighbor_ranks[r], MPI_TAG, comm, recv_requests + r);

	for (int r = 0; r < num_neighbors; ++r)
	{
		if (send_ptr[r + 1] > send_ptr[r])
		{
			int num_rows = send_ptr[r + 1] - send_ptr[r];
			int num_nonzeros = 0;
			for (int i = send_ptr[r]; i < send_ptr[r + 1]; ++i)
				num_nonzeros += Sc_LocLoc_rowptr[send_list[i] + 1] - Sc_LocLoc_rowptr[send_list[i]] + Sc_LocExt_rowptr[send_list[i] + 1] - Sc_LocExt_rowptr[send_list[i]];
			send_bufsize[r] = (num_rows + 2 * num_nonzeros) * sizeof(int);
			MPI_Send(send_bufsize + r, 1, MPI_INT, neighbor_ranks[r], MPI_TAG, comm);
		}
	}

	for (int r = 0; r < num_neighbors; ++r)
		if (recv_ptr[r + 1] > recv_ptr[r])
			MPI_Wait(recv_requests + r, &status);

	char** recv_buffer = new char* [num_neighbors];
	char** send_buffer = new char* [num_neighbors];

	for (int r = 0; r < num_neighbors; ++r)
	{
		if (recv_ptr[r + 1] > recv_ptr[r])
		{
			recv_buffer[r] = (char*)malloc(recv_bufsize[r]);
			MPI_Irecv(recv_buffer[r], recv_bufsize[r], MPI_PACKED, neighbor_ranks[r], MPI_TAG, comm, recv_requests + r);
		}
	}

	for (int r = 0; r < num_neighbors; ++r)
	{
		if (send_ptr[r + 1] > send_ptr[r])
		{
			send_buffer[r] = (char*)malloc(send_bufsize[r]);

			int position = 0;
			for (int i = send_ptr[r]; i < send_ptr[r + 1]; ++i)
			{
				int row = send_list[i];
				int length = Sc_LocLoc_rowptr[row + 1] - Sc_LocLoc_rowptr[row] + Sc_LocExt_rowptr[row + 1] - Sc_LocExt_rowptr[row];

				int* colfrom = new int[length];
				int* colind = new int[length];

				int j = 0;
				for (int k = Sc_LocLoc_rowptr[row]; k < Sc_LocLoc_rowptr[row + 1]; ++k)
				{
					colfrom[j] = rank;
					colind[j++] = Sc_LocLoc_colind[k];
				}
				for (int k = Sc_LocExt_rowptr[row]; k < Sc_LocExt_rowptr[row + 1]; ++k)
				{
					int kcol = Sc_LocExt_colind[k];
					colfrom[j] = neighbor_ranks[Sc_Ext_from[kcol]];
					colind[j++] = Sc_Ext_map[kcol];
				}

				MPI_Pack(&length, 1, MPI_INT, send_buffer[r], send_bufsize[r], &position, comm);
				MPI_Pack(colfrom, length, MPI_INT, send_buffer[r], send_bufsize[r], &position, comm);
				MPI_Pack(colind, length, MPI_INT, send_buffer[r], send_bufsize[r], &position, comm);

				delete[] colfrom;
				delete[] colind;
			}

			MPI_Isend(send_buffer[r], send_bufsize[r], MPI_PACKED, neighbor_ranks[r], MPI_TAG, comm, send_requests + r);
		}
	}

	std::map<int, int>* map_Ext = new std::map<int, int>[S2_num_neighbors];
	for (int r = 0; r < num_neighbors; ++r)
		for (int i = Sc_Ext_ptr[r], j = 0; i < Sc_Ext_ptr[r + 1]; ++i, ++j)
			map_Ext[r].insert(std::pair<int, int>(Sc_Ext_map[i], j));

	int Sc_Ext_nnz = 0;
	for (int r = 0; r < num_neighbors; ++r)
	{
		if (recv_ptr[r + 1] > recv_ptr[r])
		{
			int num_rows = recv_ptr[r + 1] - recv_ptr[r];
			int num_nonzeros = (recv_bufsize[r] - num_rows * sizeof(int)) / (2 * sizeof(int));
			Sc_Ext_nnz += num_nonzeros;
		}
	}

	int* Sc_ExtLoc_rowptr = new int[Ext_size + 1];
	int* Sc_ExtLoc_colind = new int[Sc_Ext_nnz];
	int* Sc_ExtExt_rowptr = new int[Ext_size + 1];
	int* Sc_ExtExt_colind = new int[Sc_Ext_nnz];
	int* Sc_ExtExt_colfrom = new int[Sc_Ext_nnz];

	Sc_ExtLoc_rowptr[0] = 0;
	Sc_ExtExt_rowptr[0] = 0;
	for (int r = 0; r < num_neighbors; ++r)
	{
		if (recv_ptr[r + 1] > recv_ptr[r])
		{
			MPI_Wait(recv_requests + r, &status);

			int position = 0;

			for (int i = recv_ptr[r], k = Sc_ExtLoc_rowptr[i], t = Sc_ExtExt_rowptr[i]; i < recv_ptr[r + 1]; Sc_ExtLoc_rowptr[++i] = k, Sc_ExtExt_rowptr[i] = t)
			{
				int length;
				MPI_Unpack(recv_buffer[r], recv_bufsize[r], &position, &length, 1, MPI_INT, comm);

				int* colfrom = new int[length];
				int* colind = new int[length];

				MPI_Unpack(recv_buffer[r], recv_bufsize[r], &position, colfrom, length, MPI_INT, comm);
				MPI_Unpack(recv_buffer[r], recv_bufsize[r], &position, colind, length, MPI_INT, comm);

				for (int j = 0; j < length; ++j)
				{
					if (colfrom[j] == rank)
						Sc_ExtLoc_colind[k++] = colind[j];
					else
					{
						int jfrom = r;
						if (colfrom[j] != neighbor_ranks[r])
							jfrom = map_neighbor[colfrom[j]];
						Sc_ExtExt_colfrom[t] = jfrom;
						std::pair<std::map<int, int>::iterator, bool> ret = map_Ext[jfrom].insert(std::pair<int, int>(colind[j], map_Ext[jfrom].size()));
						Sc_ExtExt_colind[t++] = ret.first->second;
					}
				}

				delete[] colfrom;
				delete[] colind;
			}

			free(recv_buffer[r]);
		}
	}

	for (int r = 0; r < num_neighbors; ++r)
	{
		if (send_ptr[r + 1] > send_ptr[r])
		{
			MPI_Wait(send_requests + r, &status);
			free(send_buffer[r]);
		}
	}

	delete[] recv_bufsize;
	delete[] recv_buffer;
	delete[] recv_requests;
	delete[] send_bufsize;
	delete[] send_buffer;
	delete[] send_requests;

	int* S2_Ext_ptr = new int[S2_num_neighbors + 1];
	S2_Ext_ptr[0] = 0;
	for (int r = 0; r < S2_num_neighbors; ++r)
		S2_Ext_ptr[r + 1] = S2_Ext_ptr[r] + map_Ext[r].size();

	int* S2_Ext_map = new int[S2_Ext_ptr[S2_num_neighbors]];
	for (int r = 0; r < S2_num_neighbors; ++r)
		for (std::map<int, int>::iterator iter = map_Ext[r].begin(); iter != map_Ext[r].end(); ++iter)
			S2_Ext_map[S2_Ext_ptr[r] + iter->second] = iter->first;

	for (int j = 0; j < Sc_ExtExt_rowptr[Ext_size]; ++j)
		Sc_ExtExt_colind[j] += S2_Ext_ptr[Sc_ExtExt_colfrom[j]];

	for (int j = 0; j < Sc_LocExt_rowptr[num_coarse]; ++j)
	{
		int r = Sc_Ext_from[Sc_LocExt_colind[j]];
		Sc_LocExt_colind[j] += S2_Ext_ptr[r] - Sc_Ext_ptr[r];
	}

	delete[] map_Ext;
	delete[] Sc_Ext_from;
	delete[] Sc_ExtExt_colfrom;

	int* S2_LocLoc_rowptr = new int[num_coarse + 1];
	int* S2_LocExt_rowptr = new int[num_coarse + 1];

	S2_LocLoc_rowptr[0] = 0;
	S2_LocExt_rowptr[0] = 0;
#ifdef USE_OPENMP
#pragma omp parallel
#endif
	{
		int* w = new int[num_coarse];
		for (int i = 0; i < num_coarse; ++i)
			w[i] = -1;

		int* v = new int[S2_Ext_ptr[S2_num_neighbors]];
		for (int i = 0; i < S2_Ext_ptr[S2_num_neighbors]; ++i)
			v[i] = -1;

#ifdef USE_OPENMP
#pragma omp for schedule(guided)
#endif
		for (int i = 0; i < n; ++i)
		{
			if (CF_map[i] != FINE)
			{
				int imap = CF_map[i];
				int cnt_loc = 0;
				int cnt_ext = 0;

				for (int k = A_LocLoc_rowptr[i]; k < A_LocLoc_rowptr[i + 1]; ++k)
				{
					if (!LocLoc_s[k]) continue;

					int kcol = A_LocLoc_colind[k];
					for (int j = Sc_LocLoc_rowptr[kcol]; j < Sc_LocLoc_rowptr[kcol + 1]; ++j)
					{
						int jcol = Sc_LocLoc_colind[j];
						if (w[jcol] != i && jcol != imap)
						{
							w[jcol] = i;
							++cnt_loc;
						}
					}

					for (int j = Sc_LocExt_rowptr[kcol]; j < Sc_LocExt_rowptr[kcol + 1]; ++j)
					{
						int jcol = Sc_LocExt_colind[j];
						if (v[jcol] != i)
						{
							v[jcol] = i;
							++cnt_ext;
						}
					}
				}

				for (int k = A_LocExt_rowptr[i]; k < A_LocExt_rowptr[i + 1]; ++k)
				{
					if (!LocExt_s[k]) continue;

					int kcol = A_LocExt_colind[k];
					for (int j = Sc_ExtLoc_rowptr[kcol]; j < Sc_ExtLoc_rowptr[kcol + 1]; ++j)
					{
						int jcol = Sc_ExtLoc_colind[j];
						if (w[jcol] != i && jcol != imap)
						{
							w[jcol] = i;
							++cnt_loc;
						}
					}

					for (int j = Sc_ExtExt_rowptr[kcol]; j < Sc_ExtExt_rowptr[kcol + 1]; ++j)
					{
						int jcol = Sc_ExtExt_colind[j];
						if (v[jcol] != i)
						{
							v[jcol] = i;
							++cnt_ext;
						}
					}
				}

				S2_LocLoc_rowptr[imap + 1] = cnt_loc;
				S2_LocExt_rowptr[imap + 1] = cnt_ext;
			}
		}

		delete[] w;
		delete[] v;
	}

	for (int i = 0; i < num_coarse; ++i)
		S2_LocLoc_rowptr[i + 1] += S2_LocLoc_rowptr[i];
	for (int i = 0; i < num_coarse; ++i)
		S2_LocExt_rowptr[i + 1] += S2_LocExt_rowptr[i];

	int* S2_LocLoc_colind = new int[S2_LocLoc_rowptr[num_coarse]];
	int* S2_LocExt_colind = new int[S2_LocExt_rowptr[num_coarse]];

#ifdef USE_OPENMP
#pragma omp parallel
#endif
	{
		int* w = new int[num_coarse];
		for (int i = 0; i < num_coarse; ++i)
			w[i] = -1;

		int* v = new int[S2_Ext_ptr[S2_num_neighbors]];
		for (int i = 0; i < S2_Ext_ptr[S2_num_neighbors]; ++i)
			v[i] = -1;

#ifdef USE_OPENMP
#pragma omp for schedule(guided)
#endif
		for (int i = 0; i < n; ++i)
		{
			if (CF_map[i] != FINE)
			{
				int imap = CF_map[i];
				int r = S2_LocLoc_rowptr[imap];
				int s = S2_LocExt_rowptr[imap];

				for (int k = A_LocLoc_rowptr[i]; k < A_LocLoc_rowptr[i + 1]; ++k)
				{
					if (!LocLoc_s[k]) continue;

					int kcol = A_LocLoc_colind[k];
					for (int j = Sc_LocLoc_rowptr[kcol]; j < Sc_LocLoc_rowptr[kcol + 1]; ++j)
					{
						int jcol = Sc_LocLoc_colind[j];
						if (w[jcol] != i && jcol != imap)
						{
							w[jcol] = i;
							S2_LocLoc_colind[r++] = jcol;
						}
					}

					for (int j = Sc_LocExt_rowptr[kcol]; j < Sc_LocExt_rowptr[kcol + 1]; ++j)
					{
						int jcol = Sc_LocExt_colind[j];
						if (v[jcol] != i)
						{
							v[jcol] = i;
							S2_LocExt_colind[s++] = jcol;
						}
					}
				}

				for (int k = A_LocExt_rowptr[i]; k < A_LocExt_rowptr[i + 1]; ++k)
				{
					if (!LocExt_s[k]) continue;

					int kcol = A_LocExt_colind[k];
					for (int j = Sc_ExtLoc_rowptr[kcol]; j < Sc_ExtLoc_rowptr[kcol + 1]; ++j)
					{
						int jcol = Sc_ExtLoc_colind[j];
						if (w[jcol] != i && jcol != imap)
						{
							w[jcol] = i;
							S2_LocLoc_colind[r++] = jcol;
						}
					}

					for (int j = Sc_ExtExt_rowptr[kcol]; j < Sc_ExtExt_rowptr[kcol + 1]; ++j)
					{
						int jcol = Sc_ExtExt_colind[j];
						if (v[jcol] != i)
						{
							v[jcol] = i;
							S2_LocExt_colind[s++] = jcol;
						}
					}
				}
			}
		}

		delete[] w;
		delete[] v;
	}

	S2.Free();

	S2.comm = comm;
	S2.LocLoc.size[0] = num_coarse;
	S2.LocLoc.size[1] = num_coarse;
	S2.LocLoc.rowptr = S2_LocLoc_rowptr;
	S2.LocLoc.colind = S2_LocLoc_colind;
	S2.LocExt.size[0] = num_coarse;
	S2.LocExt.size[1] = S2_Ext_ptr[S2_num_neighbors];
	S2.LocExt.rowptr = S2_LocExt_rowptr;
	S2.LocExt.colind = S2_LocExt_colind;
	S2.num_neighbors = S2_num_neighbors;
	S2.neighbor_ranks = S2_neighbors;
	S2.Ext_ptr = S2_Ext_ptr;
	S2.Ext_map = S2_Ext_map;
}

int ParCSR_Aggressive_HMIS_Coarsening(const ParCSR_Matrix& A, const bool* LocLoc_s, const bool* LocExt_s, int* CF_map)
{
	int n = A.LocLoc.size[0];

	int num_coarse;
	num_coarse = ParCSR_HMIS_Coarsening(A, LocLoc_s, LocExt_s, CF_map);

	ParCSR_Matrix S2;
	ParCSR_Strength2(A, LocLoc_s, LocExt_s, CF_map, S2);

	S2.SetupHalo();

	int* agg_CF_map = new int[num_coarse];

	int agg_num_coarse;
	agg_num_coarse = ParCSR_HMIS_Coarsening(S2, agg_CF_map);

	for (int i = 0, j = 0; i < n; ++i)
		if (CF_map[i] != FINE)
			CF_map[i] = agg_CF_map[j++];

	delete[] agg_CF_map;

	return agg_num_coarse;
}

int ParCSR_Aggressive_PMIS_Coarsening(const ParCSR_Matrix& A, const bool* LocLoc_s, const bool* LocExt_s, int* CF_map)
{
	int n = A.LocLoc.size[0];

	int num_coarse;
	num_coarse = ParCSR_PMIS_Coarsening(A, LocLoc_s, LocExt_s, CF_map);

	ParCSR_Matrix S2;
	ParCSR_Strength2(A, LocLoc_s, LocExt_s, CF_map, S2);

	S2.SetupHalo();

	int* agg_CF_map = new int[num_coarse];

	int agg_num_coarse;
	agg_num_coarse = ParCSR_PMIS_Coarsening(S2, agg_CF_map);

	for (int i = 0, j = 0; i < n; ++i)
		if (CF_map[i] != FINE)
			CF_map[i] = agg_CF_map[j++];

	delete[] agg_CF_map;

	return agg_num_coarse;
}

}